探索性因子分析（EFA）及其在R中实现
PCA和EFA
PCA中，通过对多元数据进行降维，将大量相关变量转化为一组很少的不相关变量，这些无关变量称为主成分（Principal Component，PC），是所有观测变量的线性组合。形成线性组合的权重都是通过最大化各主成分所解释的方差来获得，同时还要保证各主成分间不相关（彼此正交）。例如，使用PCA将30种环境变量（之间可能存在大量相关或冗余）转化为几个无关的主成分，并尽可能保留原始数据集的信息。
与PCA相比，EFA是一系列用来发现一组变量的潜在结构的方法，它通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、显式的变量间的关系。EFA的目标是通过发掘隐藏在数据下的一组较少的、更为基本的无法观测的变量，来解释一组可观测变量的相关性，这些虚拟的、无法观测的变量称为因子（Factor，F），因子被当作观测变量的结构基础或原因，而不是代表它们的线性组合。例如，使用EFA探索30种环境变量的相互关系，可用三个潜在因子（地理位置、气候环境、生物因素）来表示。
每个因子被认为可解释多个观测变量间共有的方差，因此它们常被称为公共因子。代表观测变量方差的误差（error，e）无法用因子来解释。因子和误差均无法直接测量，通过变量间的相互关系推导得到。
图注：主成分分析和因子分析模式图，展示了可观测变量（X1到X5）、主成分（PC1、PC2）、因子（F1、F2）和误差（e1到e5）。
R包psych的EFA
接下来展示R语言执行EFA的方法，以psych包的方法为例。
数据集
模拟生成了包含200个对象（行）6个变量（列）的数据集（X），接下来使用EFA探索变量间的关系。
注：真实的用于因子分析的数据中，不可能只对这么少的变量进行因子分析，此处模拟数据集仅为操作方便。
#模拟数据library(mvtnorm)
N <- 200
P <- 6
Q <- 2
Lambda <- matrix(c(0.7, -0.4, 0.8, 0, -0.2, 0.9, -0.3, 0.4, 0.3, 0.7, -0.8, 0.1), nrow = P, ncol = Q, byrow = TRUE)
set.seed(123)
Kf <- diag(Q)
mu <- c(5, 15)
FF <- rmvnorm(N, mean = mu, sigma = Kf)
E <- rmvnorm(N, mean = rep(0, P), sigma = diag(P))
X <- FF %*% t(Lambda) + E
也可首先计算6个变量间的相关性。
#计算变量间的相关矩阵，以 pearson 相关系数为例corMat <- cor(X, method = 'pearson')
head(corMat)
#相关图
library(corrplot)
corrplot(corMat, method = 'number', number.cex = 0.8, diag = FALSE, tl.cex = 0.8)
corrplot(corMat, add = TRUE, type = 'upper', method = 'pie', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
判断需提取的公共因子数量
EFA要求指定提取的公共因子数量，可在执行EFA前作个评估。
library(psych)
#确定最佳因子数量，详情 ?fa.parallel
#输入变量间的相关矩阵，并手动输入原始变量集的对象数量（n.obs）
fa.parallel(corMat, n.obs = N, fa = 'both', n.iter = 100)
#或者直接使用原始变量集
fa.parallel(X, fa = 'both', n.iter = 100)
显示提取两个因子为宜。
碎石检验的前两个特征值（蓝线三角形）都在拐角处之上，并且大于基于100次随机模拟数据矩阵的特征值均值（红色虚线）。
并且，根据Kaiser-Harris准则的特征值数大于0，也推荐选择两个因子。（注意不是根据Kaiser-Harris>1的标准，Kaiser-Harris>1是对PCA而言的，如果期望在PCA中评估选择多少个主成分合适，则可使用该标准）
提取公共因子及因子旋转
执行EFA。上述评估提示选择两个因子为宜，在该步将它作为参数输入。
#EFA 分析，详情 ?fa
#输入变量间的相关矩阵，nfactors 指定提取的因子数量，并手动输入原始变量集的对象数量（n.obs）
#rotate 设定旋转的方法，fm 设定因子化方法
fa_varimax <- fa(r = corMat, nfactors = 2, n.obs = N, rotate = 'varimax', fm = 'pa')
#或者也可直接使用原始变量集
fa_varimax <- fa(r = X, nfactors = 2, rotate = 'varimax', fm = 'pa')
fa_varimax
EFA中提取因子的方法很多，如最大似然法（ml）、主轴迭代法（pa）、加权最小二乘法（wls）、广义加权最小二乘法（gls）、最小残差法（minres）等。这里选择了主轴迭代法。
因子旋转是一系列将因子载荷阵变得更容易解释的数学方法，它们尽可能地对因子去噪。旋转方法有两种：使选择的因子保持不相关（正交旋转），和让它们变得相关（斜交旋转）。最流行的正交旋转是方差极大旋转（式中参数rotate = 'varimax'），它试图对载荷阵的列进行去噪，使每个因子只由一组有限的变量来解释（即载荷阵每列只有少数几个很大的载荷，其它都是很小的载荷）。
因子模式矩阵（标准化的回归系数矩阵，列出了因子预测变量的权重）显示，前两个因子解释了数据集的39%的方差。
对于该模拟的数据集，变量1、2在第一因子上载荷较大，变量3在第二因子上载荷较大。
使用正交旋转将人为地强制两个因子不相关。如果允许因子间存在相关，则可使用参数rotate = 'promax'。
#上述展示了正交旋转，以下是斜交旋转fa_promax <- fa(r = corMat, nfactors = 2, n.obs = N, rotate = 'promax', fm = 'pa')
fa_promax
因子模式矩阵显示，前两个因子解释了数据集的39%的方差。
因子关联矩阵显示两个因子之间存在负相关关系。在因子之间存在较大相关关系时，因子载荷阵会出现较明显的噪声（相较于正交旋转而言，因为允许潜在因子相关），虽然斜交方法更为复杂，但更符合真实数据。如果因子间的关联很低，推荐更换为正交旋转方法简化问题。
#斜交旋转中，因子结构矩阵（载荷阵）使用 F=P*Phi 获得#F 为因子载荷阵，P 为因子模式矩阵，Phi 为因子关联矩阵
fsm <- function(oblique) {
P <- unclass(oblique$loading)
F <- P %*% oblique$Phi
colnames(F) <- c('PA1', 'PA2')
F
}
fsm(fa_promax)
根据因子结构矩阵，变量1、2在第一因子上载荷较大，变量3在第二因子上载荷较大。因为这个模拟数据的两个因子之间相关程度很小，所以斜交旋转的载荷阵和正交旋转相似。
绘制正交或斜交结果的图形。
#可视化par(mfrow = c(2, 2))
factor.plot(fa_varimax, title = '正交旋转')
fa.diagram(fa_varimax, simple = TRUE, main = '正交旋转')
factor.plot(fa_promax, title = '斜交旋转')
fa.diagram(fa_promax, simple = FALSE, main = '斜交旋转')
可以据此评估哪些变量在哪些因子上的载荷更大，以确定哪些因子主要代表了哪些变量间的相互关系。
如果出现了未能被已知因子所代表的变量（例如，本模拟数据集的变量4），则表明还存在潜在的因子。如果有必要，则可以在计算过程中增添提取因子的数量。
因子得分
尽管EFA中不怎么关注因子得分，如果期望获得，需使用原始变量矩阵（而非变量间的相关矩阵）计算。
与精确计算的主成分得分不同，因子得分只是估计得到的。它的估计方法有很多种，fa()默认使用回归方法计算（计算时添加参数score='regression'）。
#因子得分，以斜交旋转结果为例fa_promax <- fa(r = X, nfactors = 2, rotate = 'promax', fm = 'pa', score = 'regression')
head(fa_promax$scores)
#得分系数（标准化的回归权重）
head(fa_promax$weight)
其它常用于EFA的R包
FactoMineR包不仅提供了PCA和EFA方法，还包含潜变量模型，并提供了分别针对数值型变量和分类变量的计算方法。FAiR包使用遗传算法估计因子分析模型。它增强了模型参数估计能力，能够处理不等式的约束条件。GPArotatio包则提供了许多因子旋转方法。nFactors包提供了用来判断因子数目的许多复杂方法。
参考资料